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ABSTRACT 

We have modified our particle-by-particle adaptation of the made-to-measure (M2M) 
method, with the aim of modelling the Galactic disc from upcoming Galactic stellar 
survey data. In our new particle-by-particle M2M algorithm, primal, the observables 
of the target system are compared with those of the model galaxy at the position of the 
target stars, i.e. particles. The mass of the model particles are adjusted to reproduce 
the observables of the target system, and the gravitational potential is automatically 
adjusted by the changing mass of the particles. This paper builds upon our previous 
work, introducing likelihood-based velocity constraints in PRIMAL. In this paper we 
apply PRIMAL to barred disc galaxies created by a iV-body simulation in a known 
dark matter potential, with no error in the observables. This paper demonstrates that 
PRIMAL can recover the radial profiles of the surface density, velocity dispersion in 
the radial and perpendicular directions, and the rotational velocity of the target discs, 
along with the apparent bar structure and pattern speed of the bar, especially when 
the reference frame is adjusted so that the bar angle of the target galaxy is aligned to 
that of the model galaxy at every timestep. 

Key words: methods: iV-body simulations — methods: numerical — galaxies: struc- 
ture — galaxies: kinematics and dynamics — The Galaxy: structure 



1 INTRODUCTION 



The made-to-m e asure 
ISver fc Tremaind (|l99e 
past few y ears, includins 
NMAGIC 



(M2M) method developed by 

has seen increasing use in the 

multiple 



involvine 



w y ears, mcludms multiple papers mvolvmg 

fe.g. Ide Lorenzi et al.|[2007l . boosl : iDas et al.|[2oTll : 

iMorganti fc Gerhard |2012| ) , and the series o f three papers 
iLong fc Mad |2010| . BoTa iLong etall l2013h bas ed on the 
M2M algorithm described in lLong fc Mad (|2O10l '). This in- 
creasing level of interest is now highlighting the potential of 
the M2M method to be used in many different ap plications, 
from probing the dark matter halo of NGC 4649 (JDas et al.l 
[201l|) to tailoring N-hoAy initial conditions (|Dehnedl2009f ). 
While previous work has prim arily focused on theor etical 
models and exte rnal gala:xies, iBissantz et al.l (|2004l ') and 
I Long et al.l ((2013) have apphed their versions of M2M to 
the Milky Way which is also our eventual goal. 

The Milky Way is k nown to be a barred spira l 
jSpergel fc Blitj|l990l : rWeTn bcre 199 HBinnev et al.lll997l ') 
but there is still disagreement over whether the bar is 
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formed of a single structure, or if it is comprised of a sep- 
arate long flat bar in addition to t h e barred bulg e (e.g . 
Martinez- Valpuesta fc GerhardI l201ll : lAthanassouJal |2012| : 



Gonzalez- Fernandez et al.ll2012l 'l. The bar angle, defined as 



an offset between the axis of the bar and the Sun- Galactic 
centre line is also still debated. For examplejpwck ct al 
( 1995) suggested a bar angle of 20° ± 10°, filcock ct J 



(|200d ) found 15 °. iRodriguez-Fernandez fc Combed (|20oi l 
found 20° - 3 5°.lMinchev et al.l l|201ll ) found 20° - 45° and 
I Green et all (|201ll 'l found 45°. There is also still discus- 
sion of the pattern speed of the bar, 0,^. Among others 
Fuxl (Il999f) cal culated a value of Q,p ~ 50 km s~^kpc~^, 
Dehned (119991) calculated flp = 53 ± 3 km s~^kpc 



Rodriguez- Fernandez fc Combed (120081') calculated fJ 



30 - 40 km s"^kpc" ^__^^^, 

60 km s~^kpc~^ and ICerhardl ( 20111 ) provided a review of 



Wanget al.l (120121 ') concluded Q.p 



previous work, concluding a likely pattern speed for the bar 
of fip ^ 50 - 60 km s^^kpc'^ 



Biss antz et al.l (|200J) were the first to apply t he orig- 
inal M2M algorithm from ISver fc Tremaind (|l99d ') to the 
Milky Way, which was an important first test of M2M on 
the Milky Way, matching observed bar kinematics in sev- 
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eral fields. The M2M algorithm has undergone significant 
improvements since then much more observational data have 
become available, leaving the door open for more up to date 
M2 M modelling of the Milky Way. 

iLong et al.l l|2013r ) ha ve taken the next st ep, applying 
their M2M algorithm from iLong fc Maol (|2010l ) to observed 
radial velocity data from the Bulge RAdial Velocity As- 
say (BRAVAJQ l|Rich et alllioOTl : iKunder et al.ll2012l ) and 
the par ticle mass distrib ution of an A'^-body barred galaxy 
model l|Shen et al.ll201Cl ) that has been built t o reproduce 
the M ilky Way disc and earlier BRAVA data. iLong et al.l 
l|2013h used the g ravitational poten tial calculated from the 
AT-body model of IShen et aD (J2010h . Then, they rotate this 
potential with a fixed pattern speed and assumed bar angle 
from the line of the Galactic centre and the observer, i.e. 
the position of the Sun. They ran many models with differ- 
ent pattern speeds and bar angles and explored the models 
which best fit their observables. They found their best model 
recove red the bar angle and pattern speed of the lShen et al.l 
l|201Cl ) A-body model, and reproduced the mean radial ve- 
locity and radial velocity dispersion of the BRAVA data very 
well. 

Although iLong et al.l (|2013l ) have taken an impor- 
tant step forwards there will be much room for improve- 
ment for applying M2M to the future observational data. 
The European Space Agency's upcoming Gaia mission 
along with ground based surveys, e.g. the Panoramic Sur- 
vey Telescope And Ra pid Response System (PanStarrs, 
e.g. iKaiser et al.l I2OICI ). the Visible and I nfra-red Sur- 
vey Telescope for Astronomy (VISTA, e.g. iMinniti et al.l 
2009|). the Larg e Synoptic Survey Telescope (LSST, e.g. 



Ivezic et al.ll2008l ). the Sloan Extension for Galactic Under - 
standing and Exploration (SEGUE, e.g. lYannv et al.|[2009^ . 
the Apache Point Ob servatory Galactic Evolutio n Experi- 
ment (AP OGEE, e.g. lAllende Prieto et all 120081 ). Subaru- 
PFS (e.g lEUis et al.l 120121'). The Radial Velocity Experi- 
ment (RAVE, e.g. ISteinmetz et al.l l2006l l. the Large sky 
Area Multi-Object fibre Spectroscopic Telescope, LAMOST 
Experiment f or Galactic Unde rstanding and Explanation 
(LEGUE, e.g. JDeng et"aLll2012l) and the Gaia -ESO public 
spectrographic survey (e.g. iGilmore et al.ll2012l ) will provide 
us with an unprecedentedly large amount of data about the 
Milky Way, which we can use as observational constraints 
for our M2M algorithm. As such we have made modifications 
to our algorithm with this goa l in mind. 

In lHunt &: Kawatal l|2013l ). hereafter Paper 1, we devel- 
oped a particle-by-particle M2M algorithm now called PRI- 
MAL (PaRtlcle- by-particle M2M ALgorithm). Because the 
Galactic stellar-survey data, such as the ones Gaia will pro- 
duce, are in the form of the position and velocity of individ- 
ual stars, PRIMAL is designed to compare the observables at 
the position of each star, i.e. not binned data as in previous 
M2M modelling. Another major difference betweenPRlMAL 
and other M2M algorithms is that the gravitational poten- 
tial is calculated via self-gravity of the model particles. The 
potential is thus altered by the changing particle masses in- 
duced by the M2M algorithm. In Paper 1 we apply primal 
to the target system of a smooth axisymmetric disc created 
by A-body simulations, and demonstrate that primal can 
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reproduce the density and velocity profiles of the target sys- 
tem well, even when starting from a disc whose scale length 
is different from the target system. 

We are encouraged by the success of Paper 1 and intend 
to build upon it. In this paper we apply primal to barred 
disc gal axies again generated by A-body simulations with 
GCD-K (|Kawata fc GibsonI l2003l : iKawata et all 120131 1. We 
introduce a new form of velo city observable constraints as 
described in lde Lorenzi et al.l (|2008). based on the likelihoo d 
function as described in iRomanowskv fc Kochanekl l|2001f ) . 
We also introduce a rotating re ference frame in a similar, 
although not identical fashion to lLong et al.l (|2013h . 

We use target systems whose information is known 
without any error. Ultimately we wish to apply primal to 
real observational data, where the information will be pro- 
vided for a limited region of the sky, with a more compli- 
cated selection and error function due to the dust extinction, 
crowding and stellar populations. However, we think that in 
the development stages it is important to test the algorithm 
against the ideal target. In this paper we demonstrate the 
successful application of primal to the barred galaxy tar- 
gets, and this is a significant step forward to modelling the 
Milky Way with M2M. 

This paper is organised as follows. Section [2.11 describes 
the M2M methodology of our original M2M algorithm as 
shown in Paper 1. Sections 12.21 and 12.31 describe the new 
improvements we have made to the algorithm, and Section 
2.4 describes the set up of our target systems. Section [3] 
shows the performance of our updated method for recreating 
the target disc systems. In Section [4] we provide a summary 
of this work. 



2 THE M2M ALGORITHM: PRIMAL 

The M2M algorithm works by calculating observable prop- 
erties from the model and the target, and then adapting 
particle masses such that the properties of the model repro- 
duce those of the target. The target can be in the form of a 
distribution function, an existing simulation, or real obser- 
vational data. The model can be a test particle simulation 
in an assumed fixed or adaptive potential, or a self-gravity 
A-body model. 



2.1 A particle- by-particle M2M 

We have presented a full description of both the origi- 
nal M2M and our particle- by-particle M2M in Paper 1. In 
this section we describe briefly the basis of our particle-by- 
particle M2M. As mentioned in Section[Tl our ultimate tar- 
get is the Milky Way, and the observables are not binned 
data, but the position and velocity of the individual stars 
which are distributed rather randomly. To maximise the 
available constraints, we evaluate the observables at the 
position of each star and compare them with the A-body 
model, i.e. in a particle-by-particle fashion. To this end pri- 
mal uses a kernel often used in Smoothed Particle Hydro- 
dynamics (SPH), W{r,h), which is a spherically symmetric 
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spline function given by 

[ 1- 

W{r,h) 



8 



Q{r/hf +&{r/hf 
2[1 - {r/h)f 




if < r//i < 1/2, 

if 1/2 ^ r/h ^ 1, 

otherwise, 

IjMonaghan fc Lattanzidll985l '). where r =| r^— r^ | and /i is a 
smoothing length described later. Note that in our particle- 
by-particle M2M algorithm the kernel, W{r,h), does not 
explicitly include the total mass, Mtot , unlike standard M2M 
algorithms, because we wish to eventually apply it to the 
Milky Way, whose mass is unknown. 

We use this kernel to calculate the density at the target 
particle locations, Tj , of both the target and the M2M model. 
For example, the density of the target at Tj is evaluated by, 



Pt,. 



= y^^mt,kW{ruj,hj), 



(2) 



where rrit^k is the mass of the target particle, r^j 
and hj is the smoothing length determined by 



Tfc- 



n 



Pt,3 



1/3 



(3) 



where ry is a parameter and we have set r; = 3. In SPH simu- 
lations, a value of r] between 2 and 3 is often used. We choose 
the higher value to maximise the smoothness. This results 
in « 113 particles being included in the smoothing length 
when the particles are distributed homogeneously in three- 
dimensional space. The solution of equation ([Sj is calculated 
iteratively until the re lative change between tw o iterations 
is smaller than 10"^ (jPrice fc Monaghanll2007l ). Similarly, 
the density at hj is calculated by 



Pj ^^m^W{^ij,hj) 



(4) 



from the model particles. The target density pt,j is calcu- 
lated only once at the beginning of the M2M simulation, 
and the model density pj is recalculated at every timestep. 
We then calculate the difference between these observables 



A, 



Pjj-t)- Pt.j 
Pt,i 



(5) 



Following the M2M algorithm, the mass of the model parti- 
cle is changed to reduce Ap- . Thus an example of the particle 
mass change equation, when using density observables is: 



dt 



mi(t) 



emi(i) 



+ p{hi 



mj: 



mi(t) 



W{rij,hj 
Pt,j 

+ 1 



■A,,p(t) 



(6) 



where rhi is the prior, and M is an arbitrary constant mass, 
which we set as M — 10^^ Mq for this paper. We introduce 
this arbitrary constant mass so that the method may be ap- 
plied to a system with unknown mass and as a result the 
parameters, such as e and /x, must be calibrated before per- 
forming the modelling. In equation (|6]) Aj,p(i) corresponds 
to a temporally smoothed version of Ap . to reduce fluctua- 
tions. This is derived by 



A.(i) 



A.(t- 



''dr. 



(7) 



with a being a small positive parameter. This Aj (t) can be 
calculated from the differential equation 

dA(i) 



dt 



= a(A- A). 



(8) 



This temporal smoothing effectively increases the number of 
particles from TV to 

r*l/2 



A^cff = iV- 



Af 



(9) 



where Ai is the length of the times tep and ti/2 = (In 2) /a is 
the half life of the ghost particles. ISver fc Tremaind (|l99d ) 
show that excessive temporal smoothing is undesirable, and 
should be limited to a ^ 2e. The regularisation term, 
pi ( In ( "'^' ' \ -|- 1 ) , in equation © is introduced to stabilise 
the mass changing process. The idea behind this term is to 
maximise the entropy defined by 



s = ~ 2_, "^» 111 



(10) 



and /i is a parameter to control the regularization. 

In this paper, the prior, rhi, is set as rhi ~ Mtat,ini/N , 
where A/tot, ini is the initial total mass of the system and A'^ 
is the number of particles. The regularization term forces 
the particle mass close to their initial value. As with Paper 
1, we write e as e = e'e" where e" is given by 

(11) 



(a^E, 



W(r,j,h,) 
Pj,t 



^pA^) 



2.2 Maximum likelihood for velocity constraints 

In Paper 1 we use velocity observables in the form of mean 
local velocity field, calculated around the target particle po- 
sitions, with the kernel described in equation ([l]). However 
as the Galactic stellar surveys will provide us velocity in- 
formation for individual particles, instead of smoothing the 
velocity, we can evaluate likelihood of the actual velocity 
of the particle. Thus we have converted the velocity section 
of our algorithm to maximise the l ikelihood of the v elocit y 
of the target particles as shown in de Lor enzi et al.l l|2008l ). 
The likelihood is calculated with the equation 



C 



J2^M^^ 



(12) 



where £, is the likelihood function for a single discrete ve- 



locity. Following iRomanowskv fc Kochanekl l|20011) . we cal 



culate the likelihood for individual velocity observables. 



at the target particle positions 
1 r fdL 



C-i ("j 



with 



Hil e-("^-"' /^"^^d^. 



(13) 



where Oj is the velocity error, which we have set as Oj — 
2.5 km s~^ for this paper, and dL/dv is a velocity distribu- 
tion for the model. Although we fix the velocity error, and do 
not discuss the effects of the errors in this paper, an advan- 
tage of the likelihood-based velocity constraints is that we 
can set individual errors for each velocity c omponent of each 
particle. Instead of the kernel chosen in Ide Lorenzi et al.l 
(l2008h we use our kernel from equation ([T]), allowing us to 
write dL/dv for target particle j, from model particle i, as 

-^Wym«5(i;j-t;0, (14) 



dv 



u 
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where 5{x) is the delta function, and 



(15) 



which is the same as equation (Q. We can express Cj in 
equation (|13[) as 






where 



and 






dnii 



(16) 



(17) 



(18) 



This leads us to the modified term in the particle mass 
change equation. Following the M2M algorithm, we max- 
imise the likelihood of equation (|12p . using 



drrii d£j 

at ami 



(19) 



where 



ACj _ d Y^ / Cj 

dm., r\m^ ^-^ \ /„ 



dm 



- E 

i 

- E 



-^ln(£,) - -^ln(;,) 
ami ami 






£, dmi 



(20) 



The particle mass change equation from the velocity based 
likelihood constraints is calculated with 



j 
and equation ([6]) becomes 

dt 



1 e-'"' 



{Vi-v,f/2a 



(21) 



mi(t) 



- CM 

+ Ew^' 

+ /i ( In 



em,{t)l A/E 



Pt,i 



A.(t) 



1 e"<"'-.3-"^,i)V2<j 



27r 



-Cr 



1 e"'"-.J"''-'*'^/^'^^.J 



r. 



27r 



]^ g-("rot,j-I'iot,i) /2CT„t J -j^ 



27r 
mi{t) 



(22) 



where Wr, tiz and Urot are the radial, vertical and rotational 
velocity components. The parameter, (, is an optional ad- 
justable parameter for changing the significance of the ve- 
locity constraints, although w e set (^ = 1 in this paper. Fol- 
lowing |deJjorenzi_et_alJ (|2008l ). we use temporally smoothed 
versions (c.f. equation[8| of £■ and /. 

Alongside this new mass change equation we have also 



altered the time when the constraints kick in from our previ- 
ous method. We found when using the likelihood-based ve- 
locity constraints the model requires a lower level of tempo- 
ral smoothing, and thus we are able to use temporal smooth- 
ing as soon as the mass change equation is enabled. Thus 
we now use the following series of stages: From t = to 
0.471 Gyr (one simulation time unit) we allow the initial 
model to experience relaxation, following a standard self- 
gravity A''-body calculation without any mass change. From 
t = 0.471 Gyr to 0.942 Gyr we used temporally smoothed 
density constraints only, and at i = 0.942 Gyr we engage 
the velocity constraints as well. This sequence is substan- 
tially shorter than the method used in Paper 1, allowing 
the solution to converge faster, and the overall simulation 
length to be halved to ~ 5 Gyr. We continue to use individ- 
ual timesteps for the particles, and only update the masses 
of particles whose position and velocity are updated within 
the individual timestep. The timestep for each particle is 
determined by 



CIZi — v^dyn 



0.5hi 



(23) 



|dVi/df|, 

with Cdyn ~ 0.2. We also retain the limit on the maxi- 
mum mass change which any particle can experience in one 
timestep. We set this limit to ten percent of that particles 
mass. 

We have again performed a parameter search to deter- 
mine differences in the likelihood-based velocity constraints, 
as we did in Paper 1. There are four important parame- 
ters, e, a, C and ^, which must be calibrated for M2M. e 
provides the balance between the speed of convergence, and 



the smoothness of the process. Note that e 



where e" 



is defined by equation (lll|l . We have chosen e' = 0.1 as an 
appropriate balance between accuracy and simulation time. 
With more computing power available to us we would con- 
sider running a lower value of e. However if e <C 0.1, it is 
possible the model might not converge as the mass change 
is too slow. The choice of a, which controls the strength 
of the temporal smoothing, should depend upon the choice 
of e (a ^ 2e). We find that our modelling is not overly 
sensitive to a as long as the condition a ^ 2e is met and 
we set a — 2.0 in this paper. We set (" = 1 as mentioned 
before. ^ controls the strength of the regularization and is 
essential in reducing the oscillation in particle masses and 
ensuring smooth convergence. We discuss the importance of 
fi in much greater detail in Paper 1. In this paper we adopt 
H — 10^. All different models presented in Section |3] use this 
same parameter set, and have not been individually tailored 
to the target or model in question. This demonstrates the 
robustness of the method. 



2.3 Rotating reference frames 

In Paper 1 it was sufficient to use a fixed reference frame as 
we were investigating smooth axisymmetric discs. However, 
if the target has some non-axisymmetric structure, such as 
a bar, the target bar angle is fixed, but the bar of the model 
rotates in the fixed reference frame. For example, if there is a 
bar, we expect the density and kinematics to be very differ- 
ent at the different azimuth angles at a fixed radius. Then, 
if the bar of the model is not aligned with the target bar, 
the observables of the model are evaluated in the different 
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dynamical states from the target observables. Hence, if the 
bar of the model keeps rotating in the fixed frame, the model 
particles receive the different constraints from the target de- 
pending on the bar angle at each timestep, and the model 
nev er settles to t he sol ution. This is discussed in Section [S] 

iLong et al.l ((2013) have proposed using a reference 
frame with a fixed bar angle, and comparing multiple simu- 
lations with different bar angles to find the best fit. This was 
trivial for their model, because they used a fixed shape of the 
bar potential and rotating with a fixed bar pattern speed, 
Qp. We however have not assumed any pattern speed prior 
to the beginning of our simulations, nor have we placed any 
explicit constraints on it. Instead we start with a smooth 
disc as the initial condition, allowing the pattern speed to 
evolve with the model galaxy due to self-gravity, and com- 
pare flp for the model and target galaxies at the end of the 
run. 

We therefore calculate the angle of the bar in our target, 
and the angle of the bar in the model at each step. Then we 
rotate the model to match the bar angle of the target for the 
purposes of calculating the observables in the same reference 
frame. It is our hope that this method will allow the pattern 
speed to be recovered along with the density and velocity 
profiles. When applying this to the Milky Way we will not 
know the exact bar angle. But here, we assume that the 
bar angle is known for our first step of modelling the bar. 
We call this reference frame change the rotating reference 
frame hereafter. In Section |3] we present a comparison of 
our method with and without this rotating reference frame, 
and also present the results from cases where we have chosen 
an incorrect bar angle. 



2.4 Target system setup 

Our simulated target galaxies consist of a pure stellar disc 
with bar and/or spiral structure and a static dark matter 
halo, set up using the method described in iGrand et al.l 
l|2012f ). The dark matte r halo density profile is taken from a 
trunc ated NFW profile (|Navarro et al.lll997l : lRodionov et al] 
I2OO9I ) and given by; 



Pdn 



3Hi 



SttG cx(1 + cx)^ 



(24) 



where Sc is the characteristic density described by 
iNavarro et al.l (|l997l ). The truncation term, e~^ , is intro- 
duced in our initial condition generator for a live halo sim- 
ulation. Although we use a static dark matter halo in this 
paper, we used the profile of equation (|24|l . Note that the 
truncation term leads to very little change in the dark matter 
density profile in the inner region focused on in this paper. 
The concentration parameter c — r2oo/rs and x — r/r^oo, 
where r^oo is the radius inside which the mean density of 
the dark matter sphere is equal to 200pcrit and given by; 



r200 = 1.63 X iO" 



Ma, 



hJoMQ 



'iiookpc, (25) 



where /iioo ~ -ffo/(fOO km s"'^ Mpc~^), and Ho is the Hub- 
ble constant set to 7f km s~^ Mpc~^. 

The stellar disc is assumed to follow an exponential sur- 
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Figure 1. Initial (red dotted), final (green dashed) and target 
(black solid) density profile (upper), radial velocity dispersion 
(upper middle), vertical velocity dispersion (lower middle) and 
rotation velocity (lower) for Model A with Target I. 



face density profile: 



Pd 



Ma ,2 

sech 



'iTTZdRl 



Zd 



-R/Rd 



(26) 



where Zd is the scale height of the disc and Rd is the scale 
length. The velocity dispersio n for each three dim ensional 
position is computed following ISpringel et al.l lJ2005h to con- 
struct a near-equilibrium condition for each the target discs. 
We have constructed four different target galaxies whose ini- 
tial conditions are listed in Table f, and the scale length 
of the target discs are initially set as Rt^d = 3 kpc. We 
run A'^-body simulations with these initial conditions, with 
10^ particles, for 2 Gyr using a tree JV-bod y code, GCD+ 
JKawata fc Gibsonll2003l : lKawata et al.l2013r ). and adopt the 
final outp ut as a target. We use the kernel softening sug- 
gested by IPrice fc MonaghanI (|2007|). Although these au- 
thors suggested adaptive softening length, we use a fixed 
softening for these simulations for simplicity. Our softening 
length e = 0.577 kpc is about three times larger than the 
equivalent Plummer softening length. We also use this soft- 
ening for the M2M modelling runs. 
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Table 1. Af-body target parameters. M200 is the mass of the halo, M^ is the mass of the disc, c is the concentration parameter, z^ is 
the scale height, a'^/a'^ is the ratio between radial and tangential velocity dispersion and Q.t,p is the pattern speed of the bar, measured 
at 2 Gyr. 



Target 


M2m (Mq) 


Ma (Mq) 


c 


Zd (kpc) 


2 / 2 


nt,p 


(km s -"^kpc 1) 


Notes 


I 


1.75 X 10^2 


3.0x101° 


20.0 


0.35 


9.0 




N/A 


Smooth Disc 


II 


2.0 X 10^2 


5.0x10^° 


9.0 


0.3 


2.0 




27.5 




III 


1.5 X 10^2 


5.0x10^0 


7.0 


0.3 


2.0 




31.7 




IV 


1.75 X 10^2 


5.0x10^° 


9.0 


0.3 


2.0 




28.9 





Table 2. M2M model results at the final timestep. Clt,p is the target pattern speed, Clp is the model pattern speed, x^ is a measure of 
accuracy of the density, Cr,z,rot arc the likelihood values for the radial, tangential and rotational velocity. 



Model 


Target 


flp (km s ^kpc 1) 


X? 


-£,/10e 


-£3/10" 


-£rot/10« 


Notes 


A 


I 




N/A 


0.123 


6.06 


3.61 


5.27 


Smooth Disc 


B 


II 




27.9 


0.254 


4.72 


3.62 


4.43 




C 


III 




30.4 


0.235 


4.96 


4.40 


5.28 




D 


IV 




28.3 


0.189 


5.14 


5.15 


5.02 




E 


II 




27.3 


0.230 


4.77 


3.65 


4.47 


Partial Data 


F 


II 




23.6 


0.276 


5.77 


3.91 


5.21 


Bar angle -30° 


G 


II 




28.0 


0.250 


4.87 


3.67 


4.55 


Bar angle -10° 


H 


II 




24.3 


0.21-0.49 


6.82-9.99 


4.59-5.48 


6.63-8.77 


No rotating frame 



As mentioned above, in this initial stage of develop- 
ment, we assume that the dark matter halo potential is 
known and there is no other external potential such as the 
bulge or stellar halo. We use the same number of particles, 
10^, and the same initial dark matter halo and disc param- 
eters for the model and target galaxies, except for the disc 
scale length: Rd = 2 kpc for the models and Rt,d = 3 kpc 
for the targets. 



3 RESULTS 

In this section we present the results from our eight models 
using PRIMAL. We will first show the results for the smooth 
featureless target disc previously explored in Paper 1. Then 
we apply primal to three different barred targets. We also 
examine how primal can reproduce the target galaxy with 
a partial data set of the observables, or an incorrect bar an- 
gle, or without the rotating reference frame using one of the 
targets. Table [2] shows which target the model is recreating, 
the bar pattern speeds, the likelihood values for radial, tan- 
gential and rotation velocity, £,, in equation (|12p and the Xp 
for the density, where 



2 

Xp 






(27) 



Note that we include only particles within 10 kpc, and Nr is 
the number of particles satisfying this criteria. Note also that 
in the likelihood case, although we seek to maximise likeli- 
hood, the values are —jC, and hence smaller values in Table 
[2] mean higher likelihood. The absolute values of jC are not 



important. We cannot compare these values between models 
for different targets. However, the relative differences in Xp 
and in jC between models for the same target observables 
are meaningful and are used in Section 3.4. 



3.1 Smooth disc 

First, we demonstrate that the newly introduced likelihood 
velocity constraints can reproduce the smooth featureless 
disc target used in Paper 1. Model A applies primal to Tar- 
get I, which was the target used in Paper 1, but using a larger 
number of particles. Note the high value of a^/a^ of Target 
I in Table 1 was used to deliberately suppress structure for- 
mation. Fig. [T] shows that the radial profiles of the density, 
the radial velocity dispersion, the vertical velocity dispersion 
and the mean rotation velocity for the target galaxy, the ini- 
tial galaxy and the final output after primal is applied. The 
figure demonstrates that primal with the likelihood-based 
velocity constraints equally or even more accurately repro- 
duces the target galaxies compared to our old version of 
the particle-by-particle M2M (Paper 1). However, a quan- 
titative comparison between the old and new version is not 
the main focus of this paper. As discussed above, we in- 
troduced the likelihood-based velocity constraints, because 
we can compare the velocity more directly and also intro- 
duce different errors for individual velocity components and 
individual particles. Therefore, the likelihood-based veloc- 
ity constraints are a necessary update, and a comparison 
with the old version is not an important issue. Note that 
the properties shown in Fig. [l] are not explicitly constrained 
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Figure 3. Same as Fig. [T] but for Model B with Target II. 



1 1 1 1 r 

Target III 




Figure 2. face-on (left) and edge-on (right) density maps of Tar- 
get II, Models B, E, F, G and H, (from top to bottom) plotted for 
comparison. The white line indicates the angle of the bar, rotated 
for comparison. The density scale is the same for all panels. 



Figure 4. Same as Fig. [2] but for Model C with Target III. 
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Figure 5. Same as Fig.[T] but for Model C with Target III. 



Figure 7. Same as Fig. [T] but for Model D with Target IV. 



1 1 1 1 r 

Target IV 




Figure 6. Same as Fig. [2] but for Model D with Target IV. 



by PRIMAL. As discussed previously in Paper 1, it is inter- 
esting to note that although our particle-by-particle ilVI2M 
uses only the first moment of the velocity components as ob- 
servables, because primal tries to reproduce the velocity of 
individual particles, the velocity distribution becomes close 
to the target and therefore the velocity dispersion can be 
reproduced as well. 



3.2 Barred disc galaxies 

In this section we present the results of Models B, C and D, 
where we apply the same parameter set for primal to model 
three different target barred galaxies. Target II is a barred 
disc galaxy showing faint spiral structure. Fig. [2] shows the 
face-on and edge-on views of Target II and the final state of 
Model B (two top panels). The final model reproduces the 
bar feature very well. The observables are only constrained 
within 10 kpc of the Galactocentric radius and hence the 
areas outside this radius are reproduced with less accuracy. 
Fig. [3] shows the radial profiles of the surface density, 
radial and tangential velocity dispersion, and the mean ro- 
tation velocity for the target and the final model compared 
to the initial model. As in Paper 1 and Model A, these ra- 
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dial profiles are not directly constrained by primal, but are 
reproduced remarkably well. Fig. [3] shows a substantial in- 
crease in radial velocity dispersion and a corresponding de- 
crease in mean rotational velocity from the initial to the 
final model. We believe that this is due to heating from the 
bar which leads to an excellent agreement with the velocity 
dispersion of the target. 

The pattern speed of the bar, fip, is also reproduced 
very well, as shown in Tables \T\ and [3 We calculate the 
pattern speed of the bar, f2p, by calculating the change in 
angle of the bar between timesteps, divided by the difference 
in time between the steps. We take Qp to be the mean value 
from the final ten steps. We found that the bar pattern speed 
of the model is Q,p — 27.9 km s~^ which is close to that of the 
target, flt,p ~ 27.5 km s~^. This is probably due to our self- 
consistent calculation of the gravitational potential, because 
once the mass distribution and kinematic properties of the 
target disc are reproduced, a bar with a similar shape and 
pattern speed to those of the target is expected to develop. 
This is certainly helped by our use of a known, fixed dark 
matter halo potential. We are pleased to see a spiral arm 
developing in the model, which looks similar to the one seen 
in the target. 

Model C applies primal to Target III, which is also a 
barred disc galaxy, but with a smaller bar than Target II 
and a boxy and peanut sh a ped b u lge (e.g . Pfennigei 



Athanassoula fc Misiriotij 20021: iDebattista et al 



1984; 



20051 : 



Bureau et al.ll2006l : ISaito et al.ll2011l ). as can be seen in the 
top panel of Fig. [l] Rather surprisingly primal reproduces 
the boxy structure of the target as shown in the bottom 
panel of Fig. [4] Fig.[5]shows the radial profiles for Model C. 
We see a slight inaccuracy in the inner 2 kpc of the radial 
velocity dispersion, and also in the rotational velocity in the 
inner 4 kpc, which corresponds roughly with the length of 
the bar. In addition, a^ is systematically higher than the 
target at all radii. As such the bar pattern speed is not as 
well reproduced as with Model B, with Q,p — 30.4 km s~^ 
compared to rit.p = 31.7 km s~^. However, we still think 
that this is a reasonably good recovery of the target, and it 
is encouraging for further developments to apply primal to 
more complicated observational data. 

Model D takes Target IV which is morphologically sim- 
ilar to Target III, with a small bar and boxy peanut feature, 
as can be seen in Fig. [S] We see a slightly larger bar in the 
model than in the target. Fig. [7] shows slight inaccuracies 
in the recovery of the radial and vertical velocity dispersion 
and mean rotational velocity in the inner 3 kpc roughly con- 
sistent with the radius of the bar. However the pattern speed 



is still recovered well with Q,p 



28.3 km s for the final 



model compared to the target of fit,p = 28.9 km s 



3.3 Working with partial data 

Even with the huge amount of data returned by Gaia and 
related stellar surveys, due to our position within the Milky 
Way's disc, we will not even come close to having a com- 
plete data set of the disc stars. Therefore it is important 
to make sure our method is still applicable when we do not 
have access to the complete picture of the disc. Our previous 
models have used all data within 10 kpc from the galactic 
centre. However Model E was performed with a simple selec- 
tion function restricting the observable volume to a 10 kpc 




2 -1 6 

Radius (kpc) 

Figure 8. Same as Fig. [T] but for Model E with Target II, per- 
formed with a partial data set. 



sphere around a point in the plane 8 kpc from the galactic 
centre, i.e. at {x, y, z) = (8, 0, 0) in Fig. [2l roughly emulating 
Gala's observable area, while ignoring effects such as extinc- 
tion and errors. This is merely the first step towards using 
primal with realistic data. 

The third panel of Fig. [2] shows the face-on and edge-on 
view of Model E, which has a similar bar to the target (top 
panel), with a hint of a spiral arm in the lower left quadrant 
matching the one visible in the target. Fig. [S] shows that an 
excellent agreement of the final model with the target radial 
profiles is still obtained with the restricted data set. This is 
an improvement on Paper 1 where we saw loss of accuracy 
when the observable field was restricted. We believe that this 
is helped by both the likelihood form of velocity observable 
and the higher resolution with which the simulations have 
been carried out. The bar pattern speed is recovered very 
well with Q,p — 27.3 km s~^ compared to the target of Q,t.p = 
27.5 km s~^. This shows the ability of primal to produce 
reasonable results when supplied with a partial data set of 
the disc particles. However we are aware that this selection 
function is crude and the next stage of our work will deal 
with more realistic selection functions and expected obser- 
vational errors. 
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Figure 9. Same as Fig. [T] but for Model F witii Target II, per- 
formed with the assumed angle of the bar offset from 30° com- 
pared to the real value. 



Figure 10. Same as Fig. [T] but for Model G with Target II, per- 
formed with the assumed angle of the bar offset from 10° com- 
pared to the real value. 



3.4 Working with an incorrect bar angle 

As mentioned in Section [Jl the bar angle of the Milky Way 
is still still debated. Ultimately, we aim to recover the dy- 
namical state of the Milky Way with primal from the future 
stellar survey data, and recovering the bar angle is also one 
of our targets. In the previous sections, we assumed that the 
bar angle of the target is known and we align the bar of the 
model galaxy to that of the target at every timestep to eval- 
uate the observables. If we do not know the bar angle of the 
target, like with the Milky Way, we could try different bar 
angles and hope that the models with the lowest x^ and/or 
the maximum likelihood values recover the bar angle of the 
target, which is the strategy taken bv lLong et al.l (|2013l ). In 
this section, we examine the effects of running primal with 
an incorrect bar angle. Models F and G are performed with 
the bar angles deliberately set to be incorrect by 30 and 10 
degrees respectively. In this section we again use all data 
within 10 kpc from the galactic centre. 

Model F has been performed while assuming that the 
bar angle is 30° less than the real angle of the target. The 
fourth panel of Fig. [2]shows a poor reproduction of the target 



bar morphology in Model F, which is significantly shorter 
than that of the target (top panel). We also see no evidence 
of the spiral structure seen in other cases. Fig. 1 101 shows the 
radial profiles for Model F. There is a discrepancy in the 
inner 4 kpc of the model compared to the target in both the 
density profile and the radial velocity dispersion. This is in 
agreement with the weaker bar shown in Fig. (2] The average 
rotational velocity is also lower across the disc. This is also 
reflected in the final pattern speed of f2p — 23.6 km s~^ 
compared to the target of Slt,p = 27.5 km s~^. However, in 
the real Milky Way case, we cannot know the correct profiles 
or the bar pattern speed in advance. On the other hand, we 
can evaluate the goodness of fit by Xp or the values of the 
likelihood, Lv In Table[21 Model F shows significantly worse 
values of Xp and £« than those of Model B which assumes 
the correct bar angle. Therefore, we should be able to tell 
easily if the bar angle is off by 30 degrees, at least in this 
simple target case. 

Model G has been performed while assuming that the 
bar angle is 10 degrees less than the bar angle of the tar- 
get. The fifth panel of Fig. [2] shows a barred disc which is 
morphologically similar to the target (top panel). The bar 
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Figure 11. Same as Fig.[T] but for Model H with Target II, which 
is performed without the rotating reference frame. 



is reproduced well whereas the spiral structure is barely vis- 
ible. Fig. [TD] shows the radial profiles for Model G, which 
again reproduces very well those of the target. The bar pat- 
tern speed is still well recovered with Q.p = 28.0 km s~^ 
compared to the target of Q,t,p = 27.5 km s~^. In Table [21 
Model G shows similar values of Xp and Cv to those of Model 
B, although the velocity likelihood values are slightly worse. 
These results may indicate that primal does not have the 
power to determine the bar angle within 10 degree accuracy, 
but can recover it with better than 30° accuracy. However, 
our ultimate target is much more complicated than this ideal 
target, and at this stage we do not explore further the ex- 
pected accuracy of recovering the correct bar angle for this 
ideal target. At least we demonstrate that with this type 
of exercise we can examine how accurately the dynamical 
model, such as primal, can recover the bar angle. In our 
future study, we will construct more realistic mock observa- 
tional data from A'^-body barred simulated discs and 'train' 
PRIMAL to recover the bar angle as accurately as possible, 
and finally evaluate the expected accuracy of our recovered 
bar angle using the comparison demonstrated in this section. 




Figure 12. Time evolution of Xp for density for Model B (black 
line) compared to Model H (red line). 



3.5 The importance of the rotating reference 
frame 

In this section we show a brief comparison between the 
resulting models with and without the rotating reference 
frame. Model H was performed under identical conditions 
to Model B, but without the reference frame corrections de- 
tailed in Section [2]3] The bottom panel of Fig. [2] shows that 
the resulting disc contains a less prominent bar, and no ev- 
idence of spiral structure in a similar fashion to Model F. 
Fig. Ill I shows the radial profiles of Model H. In the inner 4 
kpc region, the radial density and radial velocity profiles are 
lower for the model than for the target. The average rota- 
tion velocity is lower than that of the target across the whole 
disc. The pattern speed is also too low with fip = 24.3 km 
s~^ compared to the target of Q.t,p = 27.5 km s"'^. Fig. 1121 
shows a comparison of the evolution of the Xp of the density 
between Model B and Model H. The Xp in Model H expe- 
riences periodic oscillations in time with the bar rotation 
which are not seen in Model B. This lack of a smooth model 
convergence along with the poor accuracy on the recovered 
radial profiles shows the importance of having a rotating 
reference frame. 



4 SUMMARY 

We have demonstrated that our updated particle-by-particle 
M2M algorithm, primal, can recover a target disc system 
with a bar, including boxy/peanut features, in a known dark 
matter halo potential. In primal, the observables are com- 
pared with the model at the position of the target particles. 
The mass of the model particles are adjusted to reproduce 
the target observables, and the gravitational potential is cal- 
culated self-consistently from the model particle mass dis- 
tribution. We have introduced the likelihood-based velocity 
constraints to primal, which allows us to compare the ve- 
locity of the target particle more directly than the smoothed 
velocity field used in our previous algorithm. To apply this 
method to a barred disc, we evaluate at every timestep the 
density and velocity likelihood after the reference frame of 
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the model disc has been corrected, so that the bar of the 
model is always aligned with the bar of the target. Our fidu- 
cial model recovers the radial profiles of the surface density, 
the radial and vertical velocity dispersion and the mean ro- 
tation velocity of the target system very well. In addition, 
because of our self-gravity implementation of M2M, we can 
reproduce the bar morphology and pattern speed. We have 
demonstrated that primal performs well even when the ob- 
servables are restricted to within a sphere of radius 10 kpc 
around a point in the disc plane 8 kpc from the centre. 

We admit that these applications are still simplified 
cases. The excellent recovery of the properties of the tar- 
get galaxies is not surprising because we have not applied 
any observational errors or selection functions. Our ultimate 
goal is to further improve primal to be applicable to the 
future stellar survey data, including the Gaia data. While 
Gaia will return an unprecedentedly large amount of data, 
for approximately one billion stars, the accuracy of this data 
will be highly variable due to distance, extinction, location 
in the sky etc. In a forthcoming paper, we will apply pri- 
mal to more realistic mock observational data from A^-body 
simulations, taking into account the observational errors and 
selection functions. This paper also assumes a known, static, 
spherical dark matter halo potential for simplicity. In reality 
the dark matter halo remains very much unknown and yet 
has a significant effect on the dynamics of its inner galaxy. 
Thus we need to explore different dark matter halo poten- 
tials and to consider the possibility of using a model which 
includes a live halo. 

We remain optimistic that we can continue to improve 
PRIMAL, and develop a unique tool to recover the dynami- 
cal properties of the Milky Way from the future large-scale 
stellar survey data. 
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